Draft version September 20, 2010 

Preprint typeset using L£-T^X style cmulateapj v. 11/10/09 



DETECTING MASSIVE GRAVITONS USING PULSAR TIMING ARRAYS 

Kejia Lee 1,2 , Fredrick A. Jenet 3 , Richard H. Price 3 , Norbert Wex 2 , Michael Kramer 1,2 

(Dated: ) 
Draft version September 20, 2010 

ABSTRACT 

At the limit of weak static fields, general relativity becomes Newtonian gravity with a potential field 
that falls off as inverse distance, rather than a theory of Yukawa-type fields with finite range. General 
relativity also predicts that the speed of disturbances of its waves is c, the vacuum light speed, and is 
non-dispersive. 

For these reasons, the graviton, the boson for general relativity, can be considered to be massless. 
Massive gravitons, however, are features of some alternatives to general relativity. This has motivated 
experiments and observations that, so far, have been consistent with the zero mass graviton of general 
relativity, but further tests will be valuable. 

A basis for new tests may be the high sensitivity gravitational wave experiments that are now 
being performed, and the higher sensitivity experiments that are being planned. In these experiments 
it should be feasible to detect low levels of dispersion due to nonzero graviton mass. One of the 
most promising techniques for such a detection may be the pulsar timing program that is sensitive to 
nano-Hertz gravitational waves. 

Here we present some details of such a detection scheme. The pulsar timing response to a grav- 
itational wave background with the massive graviton is calculated, and the algorithm to detect the 
massive graviton is presented. We conclude that, with 90% probability, massles gravitons can be 
distinguished from gravitons heavier than 3 x 10~ 22 eV (Compton wave length A g = 4.1 x I0 12 km), 
if biweekly observation of 60 pulsars are performed for 5 years with pulsar RMS timing accuracy of 
100 ns. If 60 pulsars are observed for 10 years with the same accuracy, the detectible graviton mass 
is reduced to 5 x 10~ 23 eV (A g = 2.5 x 10 13 km); for 5-year observations of 100 or 300 pulsars, the 
sensitivity is respectively 2.5 x 10~ 22 (A g = 5.0 x 10 12 km) and 10~ 22 eV (A g = 1.2 x 10 13 km). Finally, 
a 10-year observation of 300 pulsars with 100 ns timing accuracy would probe graviton masses down 
to 3 x 10~ 23 eV (A g = 4.1 x 10 13 km). 
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1. INTRODUCTION 

Although a complete quantum version of gravitation has not yet been achieved (fSmolin 2003; Kicfcr 2006|) , in the 
weak field, linearized limit quantization of gravitation can be carried out. One can therefore ask phenomenolog icalry 
whether or not the graviton in a gravity theory is massive (IRubakovfc Tinvakovl[20M IGoldhaber fc Nietdl2^Tm. Up 
to now the most succ essful theory of gravitation is general relativity (I Will 119931; i Stairs 2003; Will 2006; lDamoul2007t 



Kramer fc Wexl [20091) . The quantization of its weak field limit (see lGuptal ()1952h and references therein) shows that 
its gravitational interaction is mediated by spin-2 massless bosons, the so called gravitons. Recently, attention has 
turned to this question due to both observational and theoretical advances . 

On the theo retical side there is the issue of the vDVZ discontinuity (jlwasakil 119701 : Ivan Dam fc Veltmanl 119701 : 
I Zakharovll 1 9 701) , the prediction that for a massive graviton, light deflection is greater than in general relativity by the 
factor 4/3, no matter how small the mass. If this were true, despite its paradoxical nature, then classical tests of starlight 
deflection would rule out massive gravitons. But recently, resolutions have been proposed to the vDVZ paradox, so 
that massive gravito ns cannot be ruled out (jVainshteinl Il972t I Visserl H~998l : iDeffavet et al.ll2"o"02t iDamour et al.l 120031 : 
iFinn fc Sutton! 120021 ) . 

On the theor etical side a l so, ad vances in higher-dimensional theories have led to graviton mass- like features (e.g., 
DGP model bylDvali et al.l (120001) 1: meanwhile, because of possible issues with stand ard gravitation from solar system 



scale ((Anderson et al.l 19981120081 ) to cosmological scale (jSanders fc McGaug h 2002), interest has been growing in al- 



terna tive gravity theories, some of which predict massive gravitons (|Dvali et al.lbodoTlAfves et al.|[2009l : lEckhardt et all 
2010). There is, therefore, considerable motivation for experimental or observational tests of whether the graviton can 
be massive. 

The literature contains many estimates for an upper limi t on the graviton mass, estimates that differ by many 
orders of mag nitude. An early lim it of 8 x 10 4 eV (lHardll973l) was based on considerations of a graviton decay to two 
photons. At about the same time IGoldhaber fc Nietol (|1974l ) used the assumption that clusters of galaxies are bound 
by more-or-less standard gravity to infer an upper limit of 2 x 10~ 29 /io eV, where ho is the Hubble constant in units 
of lOOkms^Mpc" 1 . By using the effect of graviton mass on the gen eration of gravitationa l waves by a binary, and 
the rate of binary inspiral inferred from the timing of binary pulsars, IFinn fc Suttonl (|2002|) inferred an upper limit 
of 7.6 x 10~ 20 eV. iChoudhurv et all (|2004l ) considered the effect of a graviton mass on the power spectrum of weak 
lensing; with assumptions about dark energy and other paramters, th ey estimated a upper limit o f 7 x 10~ 32 eV for the 
gravit o n mass. Rev iews about these techniques can be found, e.g. , in Gold haber fc Ni cto (2010); P article Data Group! 
(|200^ : IWill (I2O061 ). 

These upper limits are all of value, since they are based on different assumptions about the phenomenological effects 
of graviton mass. Most of the estimates are based on the effects of graviton mass on the static field of a source, 
typically assuming a Yukawa potential, though this may not be the general c ase for a particular theory of gravity, 
such as fifth force theories, or MOND-type theories (|Willl [l998: Maggiorc 2008). In view of the lack of a theory of the 
graviton, it is impor tant to have up p er lim its based on different phenomenological implications of graviton mass. 

The mass limit of IFinn fc Suttonl (|2002| ) is based on the effect of graviton mass on the generation of gravitational 
waves, not on their propagation, but the dispersion relation for propagation is also an important independent ap- 
proach to a mass limit, as has been rec ently suggested by a number of groups (jWilll 119981 : lLarson fc Hiscockl 120001 : 
iCutler et~a l. 2003; Sta vridis fc W ill 2009). Questions about this method are timely since the detection of gravitational 
waves is expected in the near future, t hanks to the progress with present ground -based laser interferometers, possi- 
ble future space-based interferometers ffloueh fc Row an 2000; [Hough et al.l I2005T ) , and pulsar timing array projects 
(jSallmen et al.lll993t iStappers et"alll2006l : lManchesterll2006l : iHobbs et al.ll2009b| ) 

The pulsar timing array is a unique technique to detect nano- Hertz gravitational waves by timing millisecond pulsars, 
which are very stable celestial clocks. It turns out that a stochastic gravita tional wave background l eaves an angula r 
dependent correlation in pulsar timing residuals for widely spaced pulsars ()Hellings fc Downall983t ILee et ahl feOOS). 
That is, the correlation C(0) between timing residual of pulsar pairs is a function of angular separation 9 between the 
pulsars. One can a nalyse the timing residual and test such a correlation between pulsar timing residuals to detect 
gravitational waves (jJenet et al.l 120051 ). We find in this paper that if the graviton mass is not zero, the form of C{9) 
is very different from that given by general relativity. Thus by measuring this graviton mass dependent correlation 
function, we can also detect the massive graviton. 

The outline of this paper is as follows. The mass of the graviton is related to the dispersion of gravitational waves 
in § HI The pulsar timing responses to a plane gravitational wave and to a stochastic gravitational wave background 
in the case of a massive graviton are calculated in § [3l The massive graviton induces effects on the shape of pulsar 
timing correlation function, which is derived in § [4J while the detectability of a massive gravitational wave background 
is studied in § [5j The algorithm for detect massive graviton using a pulsar timing array, and the sensitivity of that 
algorithm are examined in § [SI We discuss several related issues, and conclude in § 

2. GRAVITATIONAL WAVES WITH MASSIVE GRAVITONS 

We incorporate the massive gravito n into the linearized weak field theory of general relativity (jGuptal Il952t 
IArnowitt~fc Dcscr 1959: IWeinberg|ll972[ ). For linearized gravitational waves, specifying the graviton mass is equivalent 
to specifying the gravitational wave (GW) dispersion relation that follows from the special relativistic relationship 

E 2 = p 2 c 2 + m 2 c\ (1) 
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where c is the light velocity, E is energy of the particle, and p and m are the particle's momentum and rest mass 
respectively. One can derive the corresponding dispersion relation from Eq. ([1]) by replacing the momentum by p = hh g 
and the energy by E = huj g , where h is the reduced Planck constant, with k g and uj g respectively the GW wave vector 
and the angular frequency. With these replacements, the dispersion relation for a massive vacuum GW graviton 
propagating in the z direction reads 



k gK) 



W C ut 



(2) 



where e z is the unit vector in the z direction. If the gravitational wave frequency w g is less than the cut-off frequency 
uj cu t = rrigC 2 /H, then the wave vector becomes imaginary, indicating that the wave attenuates, and does not propagate. 
(The equivalent phenomena for electromagnetic waves can be found in §87 of ILandau fc Lifshitzl JT963)). 
At a space-time point (t, r), the spatial metric perturbation due to a monochromatic GW is 



h ab (t,r) =» 



E 

P=+,x 



,p p 

e ab< 



i[u!gt-r-kg(w s )] 



(3) 



where 5i indicates the real part, and where the a, b range over spacetime indicts from to 3. The summation 
is performed over the polarizations of the GW. Since we are not assuming that general relativity is the theory of 
gravitation, we could, in principle, have as many as six polarization states. For dcfinitcncss, however, and to most 
clearly show how pulsar timing probes graviton mass, we will confine ourselves in this paper to only the two standard 
polarization modes of general relativity, denoted + and x, the usual 'TT" gauge (see Appendix A [A] for the details). 
Thus, the polarization index takes on only the values P — +, x, with A p and e p standing for the amplitude and 
polarization tensors for the two transverse traceless modes. 

The polarization tensor e p is described in terms of an orthonormal 3-dimensional frame associated with the GW 
propagating direction. Let the unit vector in the direction of GW propagation be e z ; we can choose the other two 
mutually orthogonal unit vectors e x , e y to be both perpendicular to e z . In terms of these three vector, e z , e x and e y , 
the polarization tensors are given as 



— ^xa^xb *^ 



2/a 2/fc' 



*y a &xb ■ (4) 

Since the polarization tensors are purely spatial, we will have only spatial components of the metric perturbations. 
For a stochastic gravitational wave background, these metric perturbations are a superposition of monochromatic GWs 
with random phase and amplitude, and can be written as 



D_ ] v, J — OO J 



(5) 



P=+,x 



where / g = u g /2ir is the GW frequency, fi is solid angle, spatial indices i, j run from 1 to 3, and h p is the amplitude of 
the gravitational wave propagating in the direction of e z per unit solid angle, per unit frequency interval, in polarization 
state P. If the gravitational wave ba ckground is isot r opic, stationary and independently polarized, we can define the 
characteristic strain h P according to (Maggiorc 2 000t iLeee t al. 200§}), and can write 



(h p (f g ,e z )h* p '(f>,e z ')) = ^L6 PP ,S(f g - f' g )5(e z e 2 '), 



(6) 



where the * stands for the complex conjugate and () is the statistical ensemble average. The symbol 5pp> is the 
Kronecker delta for polarization states; Spp> = when P and P 1 are different, and Sppi = 1, when P and P' are the 
same. With the relationships above one can show that 



(M*)M*)> = E / 

™_ , Jo 



oo lft P,2 



(7) 



3. SINGLE PULSAR TIMING RESIDUALS INDUCED BY A MASSIVE GRAVITATIONAL WAVE BACKGROUND 

We now turn to the calculation of the pulsar timing effects due to the stochastic gravitational wave background 
prescribed by Eq. ([5]). That background is composed of monochromatic plane wave with random phase and with 
amplitude determined from some chosen spectrum. We can then first calculate the pulsar timing effect due to a 
monochromatic plane wave, and then add together all contributions to find the effects of a s tochastic background. 

In t he weak field case, the time component of the null geodesic equations for photons is ((Liang 200(J IHobbs et al.l 
l2009al) 

dui lo 2 dhij ^ „. 
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Here tu is the angular frequency for the photon, n l is the direction from the observer (Earth) to the photon source 
(pulsar), and A is an affine parameter along the photon trajectory, normalized so that in the Minkowski background 
dt/dX = u>. From this geodesic equation, the frequency shift of a pulsar timing signal induced by a monochromatic, 
plane, massive graviton GW is 



n n J 



2 (1 + (c/uj g )k g • n) 



[ft ii (t,0)-ft i3 -(t-|D|/c,D)] 



(9) 



where it is assumed that the observer is at the coordinate origin and that D is the displacement vector, in the 
background, from the observer to t he pulsar. Equation (|^) i s a m i nor gen e raliza t ion of the r elationship for zero- ma ss 
gravitons given in many references (|Estabrook fc Wahlquistl (| 1975(1 ; ISazhinl (|1978l) ; iDetweilerl (| 19791 ) ; iLee et al.l (|2008l) ) . 
It is clear that the frequency shift of the pulsar timing signal only involves the metric perturbations at the observer, 
i.e., the term /i^(i,0), and the metric perturbation at the pulsar, i.e., the term hij(t, D). The dispersion relation of 
the GWs enters, through the denominator, in the geometric factor (1 + ck g ■ n/u> g ) and well as the phase difference 
between the GW at the earth hij(t, 0) and the GW at the pulsar h{j(t— |D|/c, D). The induced pulsar timing residuals 
R(t) are given by the temporal integration of the above frequency shift at Earth, thus 



R(t) 



Aw(t) 



dr. 



(10) 



The equations Eq. ([5]), ©, and (fTTj| determine the response of pulsar timing to a monochromatic plane GW. One 
can show the pulsar timing residual R(t) induced by the stochastic GW background is 



P=+,x 



df g / dn 



(w g 



+ ck g 



-h P {f E ,e z )efAe z ) 



[1 



(11) 



This is, for example, a minor modification of Eq. (A9) given by ILee et al.l ((2(308). As we pointed out shortly before, 
the massive graviton dispersion relation enters via the term uj g + ck g • n as in Eq.([9]) and the term 1 — e _ ' ikf! ( LJs *'' D , 
which comes from the phase difference between pulsar term and earth term. 

After the nonobservable zero frequency component is removed, the auto-correlation function (R(t)R(t + r)) can be 
calculated by replacing the statistical ensemble average with the time average. 



(R(t)R(t + r)) = 



P=+,X 



,P|2 



dfg / dn 



32nf g 
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ll 



1 — cos I — (ujg + ckg • n) 



(12) 



where D = |D| is the pulsar distance. The folded timing residual power spectra SR,(/ g ) is defined to be the Fourier 
transform of the autocorrelation function of the timing residual, 



Sn(f s ) 



(R(t)R{t + T))e- 2 ^ T dT, 



for which the explicit result is 
Sr(/ s ) = 



P=+,x 



,P|2 



dil' 



167T/ g 



UJ„- 



ck„ 



1 — COS 



ckg • n) 



(13) 



(14) 



For a power-law stochastic GW background, the power spectra of the induced pulsar timing residual is (see Appendix 
for the details) 

l^(/g)| 2 



Sr(J k ) 



P=+,x 



247T 2 /. 3 



-riik 



(15) 



where the r?(/ g ) is 



and 



-4C 3 +6C+3(C 2 -l)log(^) 




if / g > /cut 
if / g < /cut 



c 



w C ut = m g c 2 /h. 



(16) 



One can see that f](f g ) is the spec tral correction f actor for massive GW backgrounds, if compared with the case of 
general relativistic GW background (|Lee et al.ll2008[ ). where rj is equal to unity. To show the frequency dependence of 
the graviton mass effect, ry(/ g ) is plotted in Fig. [TJ 
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Fig. 1. — The value of rj(f) as a function of GW frequency / in units of the cut-off frequency /cut- When the GW frequency becomes 
much larger than the cut-off frequency, r\ approaches unity, which means that the timing residuals approach those of the massless GW case. 
At the other extreme, for f g close to /cut, the power in the timing residuals is minimal. 

When f g > / cu t, a convenient polynomial approximation of ??(/) with maximal error of 1.5% is 

7?(/ g ) = -0.7764 (h?L\ + 1J48 - 1.001 (^) -0.5868^ + 1.012. (17) 

V fg J V Jg J V -tg J Jg 

A least-squares polynomial fitting technique was used to calculate the coefficients in the above equation. The RMS 
level or of the timing residual power is defined to be <r^ = J Q S^(f)df. The spectra of GW backgrounds generated by 
various astrophysical processes are usually summarized as power-law spectra with power index a, i.e. the characteristic 
strain of GWs is h c = A c (f / fo) a . For such power-law spectra, the RMS level of corresponding pulsar timing residuals 
is 

A2 f2a / f4 f3 f2 £ i \ 

[h^+h^+h^ + h^+h-a > (18) 

Jo \ Jl Jl Jl Jl Jl/ 

where the constants /3\ . . . (3$ take following values 

/?i= 3.278 x 10 _3 (a-3) _1 , 
/3 2 = 1-476 x 10" 2 (5-2a)- 1 , 
fa = 4.226 x 10 _3 (q; - 2) -1 , 
p A = 4.955 x 10" 3 (2a - 3) _1 , 

/3 5 =4.272 x 10- 3 (l-a)- x , (19) 

and the /l = Max[T _1 , / cut ] is the larger of the following two frequencies: (1) the frequency cut-off (T -1 ) due to the 
finite length time span T of observation; (2) the intrinsic frequency cut-off / cut = w C ut/ (27r) due to graviton mass. 

One can derive an u pper limit for the GW velocity using single pulsar timing data, because of the surfing effect 
(Baskaran ct al. 2008b); but it is unlikely that one can use single pulsar timing data to constrain the graviton mass 
( Baskaran et al.ir2008aD . Because of the correction factor r?(/ g ) (see Fig.[IJ, the graviton mass reduces the GW induced 
pulsar timing residuals. This prevents us from constraining the graviton mass using the amplitude of single pulsar 
timing residuals. However, as explained in the next section, the cross correlation between pulsar timing residuals from 
different directions will help us in detecting the graviton mass. 

4. THE ANGULAR DEPENDENT CORRELATION BETWEEN PULSARS 

A stochastic GW background leaves a correlation between timing residuals of pulsars pairs HHellings fc Downslll98l 
iLee et aLll2008D . Such a correlation, C(0), depends on the angular distance 9 between two pulsars. It turns out that 
the graviton mass changes the shape of this correlation function. One can therefore detect a massive graviton by 
examining the shapes of pulsar timing correlation functions. 

As shown in Appendix [B] the pulsar timing cross-correlation function for a massive GW background depends on 
the graviton mass, specific power spectra of GW background, and observation schedule. In this way, an analytical 
expression for the cross-correlation function would not be possible. We use Monte-Carlo simulations in this paper 
to determine the shape of the correlation function for GW backgrounds with a power-law spectra. In the Monte- 
Carlo simulations for C(8), we randomly choose pulsars from an isotropic distribution over sky positions. We then 
hold constant these pulsar positions and calculate the angular separation 9 between every pair of pulsars. Next, to 
simulate the power-law GW background, we generate 10 4 monochromatic waves, choosing random phase, and choosing 



(> 



the amplitude from the power-law h c = A c (f/fo) a , where we take a = -2/ 3 (jPhinnevi |200H Uaffe fc Backerl [20031: 
iWvithe fc Loe b 2003; Enoki et al.l 120041: iSesana et all 12004 1 Wen et al.ll2009fh The timing residuals are calculated 
using Eq. (J9j) and (|T0| . Then, the cross-correlation function, C(9), between pulsar pairs is calculated. We repeat such 
processes and average over the angular dependent correlation function C(9) until the change in C(9) is less t han 0.1%. 
The a veraged correlation function is then smoothed by fitting an eighth order Legendre polynomial (see I Lee et aLI 
(2008) for the details). This smoothed C{9) is the correlation function we need. The C(9) are plotted for various 
parameters in Fig. [5J We also check the results by choosing different sets of pulsars to make sure the C{9) is not 
sensitive to the details of the random pulsar samples. As the C{9) is the statistical expectation of the correlation 
function, we call it the theoretical correlation function in contrast with the observed correlation function defined in the 
next section. 

As one may expect, the massive graviton has stronger effects for data from long observing periods than from short 
periods. One can see this by comparing the 5-year and 10-year correlation functions given in Fig. [2j where curves 
corresponding to the same range of graviton mass show considerably greater deviations from the massless case in the 
10-year correlation function than in the 5-year one. 

5 yr 1 yr 




30 60 90 120 150 30 60 90 120 150 

e e 



Fig. 2. — The atlas for cross-correlation functions C(9). The label of each curve indicates the corresponding graviton mass in units of 
electron-volts (eV). The left panel shows the correlation functions for a 5 year bi-weekly observation. The right panel shows correlation 
functions for 10-years of bi-weekly observations. We take a = —2/3 for these results. These correlation are normalized such that the 
C(0) = 0.5 for two different pulsars. 



5. AN ESTIMATE OF THE DETECTABILITY OF A GRAVITATIONAL WAVE BACKGROUND WITH A MASSIVE GRAVITON 

As we have explained, the massive graviton reduces the pulsar timing response to GWs through the correction 
factor r](f g ). In this way, a non-zero graviton mass reduces the pulsar timing array sensitivity for detecting a GW 
background. It is interesting to know how the sensitivity changes, if the graviton is massive. 

The stochastic GW background is detected by comparing the measured cross-correlation function c(9i) with the 
theoretical correlation C{9) calculated in previous section. Here, the measured cross-correlation function c(9i) is 
defined by 

' ^N-l/r, ,. \ n / . \ \ 9 



where the R a (ti) and Rb(ti) are the timing residuals of pulsar 'a' and 'b' at time ti and N is the number of observations. 

The Rb{ti) = Y^q 1 Rb{ti)/N and R a {U) = Ezlo 1 Ra{U)/N. The 9 m is the angle between the direction pointing to 
pulsar 'a' and the direction pointing to pulsar 'b'. Given N p pulsars, the index m runs from 1 to the number of pulsar 
pairs M = ( N„ — l)N n /2 beca use the autocorrelations are not used. 
Following Uenet et all (|2005| ). we define 

E" =1 (C{0 m )-C)(c(0 m )-c) 



E™=i (C(0 m ) - C) 2 E'Li (c(9 m ) - cf 

where C — ^ j C(9 m )/M and c = Em=i C (@m)/M- Then the statistic S, describing the significance of the detection, 
is S = vM p. In particular, when there is no GW present, the c(9 m ) will be Gaussian-like white noise, the probability 
of getting a detection significance larger than S is about erfc(S/\/2)/2 ()Jenet et a l. 2005). 

Our aim is to determine the ability of a given pulsar timing array configuration to detect a GW background. To do 
this, we calculate the expected value for the detection significance S by using a second set of Monte-Carlo simulations. 
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These second Monte-Carlo simulations are similar to the first ones, but instead of calculating the average value for 
C(9), we inject white noise for each pulsar, to represent the intrinsic pulsar noise and instrumental noise, and we 
calculate the expected value of S. We summarize the steps here. 

1. Generate a large number of GW source (10 4 ) to simulate the required GW background. 

2. Calculate the timing residual for each pulsar as described above, and add white Gaussian noise. 

3. Calculate the measured correlation c(6 m ) using Eq. pO)) . and calculate the detection significance S using Eq. ([21]) . 

4. Repeat steps 1,2,3 and average over the detection significance S. The converged S is the value needed to estimate 
the detection significance. 

The results for the expectation value of S, as a function of GW amplitude A c for various pulsar timing array 
configurations, are presented in Fig. [3] We have also compared simulations from several different pulsar samples with 
the same number of pulsars to make sure such S is not sensitive to the detailed configuration of the pulsar samples. 

Two features of the curves in Fig. [3] are worth noting. First, the minimal detection amplitude of a GW background 
becomes larger, when a massive graviton is present, i.e., the leading edge of the S-A c curve shifts rightwards as m g is 
made larger. This tells us that in order to detect a massive GW background one needs a stronger GW background 
signal or a smaller pulsar intrinsic noise than in the case of a massless GW background. As previously noted, this effect 
is mainly due to the reduction of the pulsar timing response and the reduction of the gravitational wave amplitude at 
lower frequencies. Fig. |3] also tells us when we can neglect the effect of a massive graviton. It is clear from Fig. [3] that 
if TO g < 2 x 10 -23 eV for a 5-year observation, the minimal detection amplitude is not reduced by more than 5%. For 
10 years of observation a 5% reduction corresponds to m g = 10 -23 eV. 

The second noteworthy feature of the S-A c curves in Fig. [3] is that of the saturation level of detection significance. 
Due to the pulsar distance term of Eq. (JTTJ) (the term involving the D), the detection significance achieves a saturation 
level when the GW induced timing residuals are much stronger than the intrinsic pulsar timing noise (jJenet et al.l 
2005). From Fig. El we note that the saturation level of detection significance is large, when the graviton is massive, 
i.e., the plateau at the right part of the S-A r curve be comes higher for a more massive graviton. This is rather similar 
to the whitening filter discussed by (jJenet et a l. 2005). The graviton mass introduces a low frequency spectral cut-off, 
which is equivalent to applying a whitening filter to the timing signals. In this way, the saturation level of S starts to 
grow for a massive GW background. Because the cut-off in the frequency domain coherently removes low frequency 
G W components, the envelope of the these S-A c curves is similar to the curves with the whitening filter as described 
bv lJenet et all (|2005| ). 
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Fig. 3. — The expected GW background detection significance using a pulsar timing array with 20 pulsars, observed for 5 years, with 
100 ns timing noise. The graviton mass, in units of electron volts, is labeled above each curve. The x-axis is the amplitude for the 
characteristic strain of the GW background (/o = lyr - 1 , a = —2/3), while the y-axis is the expected detection significance S. 



6. ESTIMATE OF THE DETECTABILITY OF THE MASSIVE GRAVITON 

In §[5l we have discussed how the detection significance for the stochastic GW background is affected by the graviton 
mass. Here we formulate the qustion as a detection problem rather than a parameter estimation problem. In this 
section, we will construct a detector, which accepts pulsar timing array data and then determines if the gravitational 
wave background is massless or massive. Then we simulate pulsar timing data, passing through the detector, to 
determine the quality of such detector. With the detector quality, we then discusse related technical requirements (the 
number of pulsars, the pulsar noise level, and so on) for performing such detections. 

The question, whether the graviton mass is zero or not, is best formulated as a statistical hypothesis test composed 
of two hypotheses F and H given by 

J H: The graviton is massless; 
1 F: The graviton is massive. 




The detection algorithm determines which hypothesis is accepted. We define the detection rate Pd as the probability 
that the detection algorithm gives statement F, when the graviton is massive, i.e., Pd = P(F\ m g > 0); and we define 
the false alarm rate Pf as the probability of getting statement F when the graviton is massless, i.e., Pf = P(F\m g = 0). 
The quality of the detection algorithm is evaluated by calcu lating the relation betwee n the detection rate Pd and the 
false alarm rate Pf. Here we take the standard approach (jDiFranco fc Rubml Il968f ) that one fixes Pf to a certain 
threshold level Pth and calculates the corresponding Pj. Throughout this paper we fix P t y, = .1%. If a detector 
maximizes the Pd for a prescribed value of Pf , we say that such detector is optimal (Kassam 1988). 

As we have shown in §21 the graviton mass changes the s hape of the correlation function. The best way to detect 
such a difference is to use following statistics (jKassa 

= Em=l (Cm(flm) ~ KM ~ t) 

7 vSEi ( c ^ dm) - ^5 (c(0m) 1 ^ v^S ( c °^ - ^) 2 e"U - C) 2 ' 

where Cq(9) is the correlation function for a massless GW background, and C m (6) is the correlation function for a 
massive GW background, which maximizes the detection significance S among possible t heoretical corre lation functions 
for all possible values of m g . One can show that the optimal statistical decision rule is (|Kassa 

J Choose H: if 7 < 7th, 
1 Choose F: if 7 < 7th, 

where the 7th is the threshold of statistical decision, which is determined to guarantee the false alarm rate constraint 
Pf<Ph = 0.1%. 

We determine the threshold 7 t h by a third set of Monte-Carlo simulations. These simulations take following steps. 
First, a massless GW background is generated, then we search for the matched value of m g , such that the corresponding 
correlation function C m (9) maximizes the detection significance S. Then we use this C m ((9) to calculate the statistic 
7. We repeat this, recording the values of 7, to establish the statistical distribution of 7 values. We take the value of 
the threshold 7th to be that for which there is less than a 0.1% chance of getting 7 > 7th for the case of m g = 0. 

After we determined the 7th, a fourth Monte-Carlo simulation is used to calculated Pi- This simulation is very 
similar to the previous one, except that a massive GW background is generated. We repeat the simulation 10 3 times 
and take as Pj the probability of getting 7 > 7th- 

We summarize the results of the simulation in Fig. HI which are gray-scale contour plots for the detection rate Pd for 
different scenarios of using 60, 100, and 300 pulsars respectively. The corresponding parameters are given in the legend 
of each panel. Intuitively, the necessary conditions for a positive detection of a graviton mass should be first, that the 
GW is strong enough for the GW to be detected, and second, that the graviton mass is large enough to change the 
shape of correlation function. This intuition is confirmed by our simulations, which show that the high detection rate 
concentrates in the upper right corner of each panel, where both graviton mass and GW amplitude are large enough. 

From Fig. |4] we can also see that we need at least 60 pulsars to be able to tell the difference between a massive GW 
background and a massless one. For 5-year observations of 100 pulsars we can start to detect a graviton heavier than 
2.5 x 10~ 22 eV and we can achieve a limit of m g = 10 -22 eV by using 5-year observations of 300 pulsars. We can 
achieve levels of 10 -22 eV and 5 x 10~ 23 eV in 10-year observations using 100 and 300 puslars respectively. 

We also note that there is a positive correlation between the minimal detectable graviton mass and the GW back- 
ground amplitude, as shown by the leftwards lead edge of the contours, in other words, we need a larger GW amplitude 
such that the GW background can be detected, if gravitons are more massive. As discussed above, this correlation is 
due to the reduction of pulsar timing residuals due to massive graviton. 

7. CONCLUSION AND DISCUSSION 

Current efforts to detect gravitational waves using radio pulsar observations make use of correlated arrival time 
fluctuations induced by the presence of a stochastic GW background. Einstein's general theory of relativity makes 
a specific prediction for the angular dependence of the correlation. The power spectrum of the induced fluctuations 
is given by the power spectrum of the gravitational wave background, determined by the physical processes of the 
generation process, divided by the square of the GW frequency. In this paper, we showed that the form of the expected 
correlation function, as well as the power spectrum of the induced timing fluctuations, will be significantly altered if 
the graviton had a non-zero mass. Only the transverse traceless GW modes were considered in this analysis. 

A non-zero graviton mass introduces a cut-off frequency below which no GWs will propagate. If we consider a fixed 
amount of time over which the correlation will be measured, five years for example, a small increase in the graviton 
mass will actually increase our ability to detect a strong background. Here, we are assuming that the background is 
being generated by an ensemble of supermassive black hole binaries. This effect is due to the fact that the presence 
of the graviton mass will act to flatten out the spectrum of the induced time residuals (see Eq. (fTB"]) and Fig. [IJ, thus 
making the strong background easier to detect. As the mass is increased, the cut-off frequency will increase. Once 
this cut-off frequency is above the inverse of the observing time, the sensitivity starts to fall of dramatically. This is 
due to the fact that pulsar timing is most sensitive to the lowest observable frequencies and the presence of a massive 
graviton removes all gravitational wave power at frequencies below the cut-off frequency. 

A non-zero graviton mass will also change the shape of the angular dependent correlation function. In order to 
understand this, consider two GW waves traveling near the direction of the line of sight between the Earth and a given 



Z2=l (Co(flm)-Co) (C(fl„ 
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Fig. 4. — Contours of 50%, 70 % and 90% detection rates for graviton mass. The white areas are the parameter space for less than 50% 
detection rate; the light gray regions are the parameter space with more than 90% chance of detecting graviton mass. The false alarm 
rate is fixed at 0.1% for all calculations. The time span between two successive observation is two weeks. The horizontal axis is the base 
10 logarithm of the characteristic strain for the GW background. The vertical axis is the logarithm of graviton mass in unit of eV. The 
number of pulsars in the timing array, the total time span of observation, and the level of intrinsic pulsar timing noise are given in the 
legend of each panel. For all the results, we use a = —2/3. 

pulsar. One GW is traveling towards the pulsar, and the other is traveling towards the Earth. The induced timing 
fluctuations will be very different for each of these waves. The GW traveling towards the Earth will induce higher 
amplitude fluctuations then its counter part traveling towards the pulsar. Now, consider the same case but with a 
massive graviton. As the GW frequency approaches the cut-off frequency, the GW wavelength get arbitrarily large. 
Assume that the GW wavelength is much larger than the Earth-pulsar distance. In this case, both GWs look exactly 
the same as far as the Earth-pulsar detector is concerned. Hence, both waves induce the same timing fluctuations. 
This symmetry between the two GW directions implies that the timing fluctuations from two pulsars near each other 
on the sky or 180° apart will have the same amount of correlation, unlike the massless graviton case. Thus, as the 
graviton mass increases, the correlation curve will become more symmetric about 90°. 

Since the graviton mass changes the expected correlation function, we can detect the massive graviton by measuring 
the shape of correlation function. The optimal detection algorithm differentiating between a massive and a massless 
GW background was constructed in this paper. Using this algorithm, we found that at least 60 pulsars are required 
to discriminate between a GW background made up of massive a nd massless gravitons. Note that, like the case of 
discriminating various polarization modes of GWs (|Lee et al.ll2008[ ). the minimum number of pulsars is insensitive to 
the intrinsic timing noise level or the total observing time. 

Here, we are mainly focusing on answering following question: What is the threshold of graviton mass, such that 
the detector, we constructed in the paper, will find out that the GW background is composed of massive gravitons 
rather than masslessgravitons? We answered this statistical detection question in this paper. The related parameter 
estimation question □ will be investigated in the future works. 

In this paper, we treat the graviton mass in the sense of a GW dispersion relation. A better graviton mass upper limit 
(m R < 2x 10~ 26 e V) can be achieved by observing the GW dispersion using Laser Interferometer Space Antenna (LISA) 
(IStavridis fc WiHll2009h . Besides GW disper sion based techniques, there are other good upper limits (2 x 10 29 /in eV) 



4 How well can we measure the graviton mass or put uplimit on the graviton mass? 
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from galaxy cluster observations (jGoldhaber & Nicto 1974). However since some gravity theories dMaggiorel 12008ft 
contain mass terms for h a b, but maintain the scalar sectors, the m g upper limits from GW dispersion observations are 
independent of the upper limits from Yukawa potential experiments such as the Solar System and the galaxy cluster 
observations. 

In summary, for the task of detecting the massive graviton using a pulsar timing array, there is one critical re- 
quirem ent: a large sample of stable pulsars. Thus the on-go ing and upcoming pr ojects like the Parkes Pulsar T iming 
Array (jHobbs et al.ll2009bO . Europ ean Pulsar Timing Ar ray (|Stappers et al.ll2006[ ). NANOGrav (j Jenet et al.ll2009D . the 
Large European Array for Pulsars dStappers et al.ll2009l) . the Five- hundred- meter Aperture Spherical Radio Telescope 
(|Nan et al.ll2006l ; iSmits et"al"1l2009D and the Square Kilometer Array (SKA) will offer unique opportunities to find, to 
time these pulsars; and to detect the GW background and measure its properties. 

8. ACKNOWLEDGMENT 

K. J. Lee gratefully acknowledge support from ERC Advanced Grant "LEAP", Grant Agreement Number 227947 
(PI M. Kramer), F. A. Jenet and R. H. Price gratefully acknowledge support from the NSF under grants AST0545837 
and PHY0554367, and support from the UTB Center for Gravitational Wave Astronomy. We also thank M. Pshirkov 
and N. Straumann for very useful comments. 

REFERENCES 



Alves, M. E. S., Miranda, O. D., & de Araujo, J. C. N. 2010, 

CQG,27, 145010 
Anderson, J. D., Campbell, J. K., Ekelund, J. E., Ellis, J., k 

Jordan, J. F. 2008, Physical Review Letters, 100, 091102 
Anderson, J. D., Laing, P. A., Lau, E. L., Liu, A. S., Nieto, M. M., 

k Turyshev, S. G. 1998, Physical Review Letters, 81, 2858 
Arnowitt, R., k Deser, S. 1959, Physical Review, 113, 745 
Baskaran, D., Polnarev, A. C, Pshirkov, M. S., k Postnov, K. A. 

2008a, Phys. Rev. D, 78, 089901 
— . 2008b, Phys. Rev. D, 78, 044018 

Choudhury, S. R., Joshi, G. C, Mahajan, S., k McKellar, 

B. H. J. 2004, Astroparticle Physics, 21, 559 
Cutler, C, Hiscock, W. A., k Larson, S. L. 2003, Phys. Rev. D, 

67, 024015 

Damour, T. 2007, A Century from Einstein Relativity: Probing 
Gravity Theories in Binary Systems: Binary Systems as 
Test-beds of Gravity Theories, ed. M. Colpi (Springer) 

Damour, T., Kogan, I. I., k Papazoglou, A. 2003, Phys. Rev. D, 
67, 064009 

Deffayet, C, Dvali, G., Gabadadze, G., k Vainshtein, A. 2002, 

Phys. Rev. D, 65, 044026 
Detweiler, S. 1979, ApJ, 234, 1100 

DiFranco, J. V., k Rubin, W. L. 1968, Radar Detection, ed. 

W. L. Everitt, Prentice-Hall Electrical Engineering Series 

(Prentice-Hall INC) 
Dvali, G., Gabadadze, G., k Porrati, M. 2000, Physics Letters B, 

485, 208 

Eckhardt, D. H., Pestana, J. L. G., & Fischbach, E. 2010, New 

Astronomy, 15, 175 
Enoki, M., Inoue, K. T., Nagashima, M., k Sugiyama, N. 2004, 

ApJ, 615, 19 

Estabrook, F. B., & Wahlquist, H. D. 1975, General Relativity 

and Gravitation, 6, 439 
Finn, L. S., & Sutton, P. J. 2002, Phys. Rev. D, 65, 044022 
Goldhaber, A. S., k Nieto, M. M. 1974, Phys. Rev. D, 9, 1119 
— . 2010, Reviews of Modern Physics, 82, 939 

Gupta, S. N. 1952, Proceedings of the Physical Society A, 65, 161 
Hare, M. G. 1973, Canadian Journal of Physics, 51, 431 
Hellings, R. W., & Downs, G. S. 1983, ApJ, 265, L39 
Hobbs, G., Jenet, F., Lee, K. J., Verbiest, J. P. W., Yardley, D., 
Manchester, R., Lommen, A., Coles, W., Edwards, R., k 
Shettigara, C. 2009a, MNRAS, 394, 1945 
Hobbs, G. B., Bailes, M., Bhat, N. D. R., Burke-Spolaor, S., 
Champion, D. J., Coles, W., Hotan, A., Jenet, F., 
Kcdziora-Chudczer, L., Khoo, J., Lee, K. J., Lommen, A., 
Manchester, R. N., Reynolds, J., Sarkissian, J., van Straten, 
W., To, S., Verbiest, J. P. W., Yardley, D., k You, X. P. 2009b, 
Publications of the Astronomical Society of Australia, 26, 103 
Hough, J., k Rowan, S. 2000, Living Reviews in Relativity, 3, 3 
Hough, J., Rowan, S., k Sathyaprakash, B. S. 2005, Journal of 

Physics B Atomic Molecular Physics, 38, 497 
Iwasaki, Y. 1970, Phys. Rev. D, 2, 2255 
Jaffe, A. H., & Backer, D. C. 2003, ApJ, 583, 616 



Jenet, F., Finn, L. S., Lazio, J., Lommen, A., McLaughlin, M., 
Stairs, I., Stinebring, D., Verbiest, J., Archibald, A., 
Arzoumanian, Z., Backer, D., Cordes, J., Demorest, P., 
Ferdman, R., Freire, P., Gonzalez, M., Kaspi, V., Kondratiev, 
V., Lorimer, D., Lynch, R., Nice, D., Ransom, S., Shannon, R., 
k Siemens, X. 2009, Preprint, arxiv: 09091058 

Jenet, F. A., Hobbs, G. B., Lee, K. J., k Manchester, R. N. 2005, 
ApJ, 625, L123 

Kassam, S. A. 1988, Signal Detection in Non-Gaussian Noise 

(New York: Springer Verlag) 
Kiefer, C. 2006, Annalen der Physik, 15, 129 

Kramer, M., k Wex, N. 2009, Classical and Quantum Gravity, 26, 
073001 

Landau, L. D., k Lifshitz, E. M. 1960, Electrodynamics of 

continuous media (Oxford: Pergamon Press) 
Larson, S. L., k Hiscock, W. A. 2000, Phys. Rev. D, 61, 104008 
Lee, K. J., Jenet, F. A., & Price, R. H. 2008, ApJ, 685, 1304 
Liang, C. B. 2000, Introduction of differential geometry and 

general relativity (in Chinese), Vol. VI (Beijing: Publisher of 

Beijing Normal University) 
Maggiore, M. 2000, Phys. Rep., 331, 283 
— . 2008, Gravitational waves, Vol. VI (Oxford: Oxford 

University Press) 
Manchester, R. N. 2006, Chinese Journal of Astronomy and 

Astrophysics Supplement, 6, 020000 
Nan, R., Wang, Q., Zhu, L., Zhu, W., Jin, C, k Gan, H. 2006, 

Chinese Journal of Astronomy and Astrophysics Supplement, 6, 

020000 

Particle Data Group. 2008, Physics Letters B, 667, 1 
Phinney, E. S. 2001, preprint, arxiv: 0108028 

Rubakov, V. A., k Tinyakov, P. G. 2008, Physics Uspekhi, 51, 759 
Sallmen, S., Backer, D. C, k Foster, R. 1993, in Bulletin of the 

American Astronomical Society, Vol. 25, Bulletin of the 

American Astronomical Society, 853 
Sanders, R. H., k McGaugh, S. S. 2002, ARA&A, 40, 263 
Sazhin, M. V. 1978, Soviet Astronomy, 22, 36 
Sesana, A., Haardt, F., Madau, P., k Volonteri, M. 2004, ApJ, 

611, 623 

Smits, R., Lorimer, D. R., Kramer, M., Manchester, R., Stappers, 

B., Jin, C. J., Nan, R. P., k Li, D. 2 009, A&A, 505, 919 
Smolin, L. 2003 , |arXiv:hep-th/0303185| 
Stairs, I. H. 2003, Living Reviews in Relativity, 6, 5 
Stappers, B., Vlemmings, W., k Kramer, M. 2009, Proceedings of 

Science, EXPReS09, 20 
Stappers, B. W., Kramer, M., Lyne, A. G., D'Amico, N., k 

Jessner, A. 2006, Chinese Journal of Astronomy and 

Astrophysics Supplement, 6, 020000 
Stavridis, A., k Will, C. M. 2009, Phys. Rev. D, 80, 044002 
Vainshtein, A. I. 1972, Physics Letters B, 39, 393 
van Dam, H., k Veltman, M. 1970, Nuclear Physics B, 22, 397 
Visser, M. 1998, General Relativity and Gravitation, 30, 1717 
Weinberg, S. 1972, Gravitation and Cosmology: Principles and 

Applications of the General Theory of Relativity (John Wiley 

k Sons) 



11 

Wen, Z. L., Liu, F. S., & Han, J. L. 2009, ApJ, 692, 511 Wyithe, J. S. B., & Loeb, A. 2003, ApJ, 590, 691 

Will, C. 2006, Living Reviews in Relativity, 9, 3 Zakharov, V. I. 1970, JETP, 12,312 

Will, C. M. 1993, Theory and Experiment in Gravitational 

Physics (UK: Cambridge University Press) 
— . 1998, Phys. Rev. D, 57, 2061 



12 



APPENDIX 
POWER SPECTRA OF TIMING RESIDUAL 
The power spectra of timing residuals is calculated by integrating Eq. (fTl"]) . Here we give the details of the calculation. 



Pulsar 1 




Fig. 5. — The geometric configuration of the coordinates and unit vectors used here. The X, Y and Z are the coordinate unit vectors, 
e z is the propagation direction of the GW, and ni and n2 are unit vectors pointing to the pulsars. 



In Fig. O the components of e z in the X a = (X, Y,Z) frame can be seen to be (sin 9 g cos <f> g , sui9 g sin</> 9 , cos# 9 ), 
where the 9 g and <f> g are respectively the polar angle and azimuthal angle of the GW propagation vector. To proceed, 
we need the components of the polarization tensors in the X a frame. The transformation from the components eo, c d 
given in the GW frame = (e a , e H , e z ) of Eq. (j4j is made with e^ b = T ca Tdb £q c d> where T ca = e c • X Q has components 



(cos 9 g cos 4> g cos ipg — sin (j> g sin ipg cos 6 g cos ip g sin <p g + cos <p g sin ip g — cos ip g sin g \ 
— cos ip g sin <f) g — cos 9 g cos (p g sin ip g cos 4> g cos if) g — cos 9 g sin cp g sin ip g sin 9 g sin ip g 
cos g sin g sin 9 g sin g cos 9 g ) 

Since the GW background is isotropic, with no loss of generality one can choose n = {0, 0, 1} so that 

w g + cn • k g = (1 + £cos6> g ) Wg) 

where 



c 



T 



'l 



J cut 

72" 

Jg. 



(Al) 

(A2) 
(A3) 



for / g > /cut • The term efjh 1 ^ now simplifies to 

e+n'n 3 = sin 2 g cos(2^ g ) 
eyn*n j = - sin 2 6> g sin(2^ g ) . 
After the ensemble average over the polarization angle ip g , the integration of Eq. (fT4|) gives 



E 

P=+,x 



1287T 3 /| 
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sin 



(l + Ccos6» g ) 2 



1 — cos 



Duj a 



(1 + Ccosflg) 



(A4) 



(A5) 



For practical pulsar timing array, we have $0 = Du g /c 3> 1, i.e. , the GW wave length is much smaller than the 
pulsar-earth distance, so one gets 



1 (1 - u 2 ) 2 4 
' [l-cos($o(l + CM))]^=- 



-4C 3 + 6C + 3(C 2 -1) log(g) 



U (1 + C/i) 2 L v uv ' 3C 5 

Finally, the power spectra of the pulsar timing residuals then becomes 

l^(/g)| 2 



o (% 2 ) 



P=+,x 



24^2/1 



(A6) 



(A7) 
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For / g > / cut , r) is given by 

-4C 3 + 6C + 3 (C 2 - 1) log 
*K/g) = ^ ' (A8) 

with C is given by Eq. (|A3[) . For / g < / cut the gravitational waves cannot propagate, and we take 77 = 0, 

ANGULAR DEPENDENT CORRELATION FUNCTION FOR PULSAR TIMING RESIDUALS 
The correlation function C(6) between pulsar 'a' and '6' defined in this paper is 

where 9 is the angle between the direction to pulsar 'a' and the direction to pulsar The R a (t) and Rb(t) are the 
GW induced timing residuals for puls ar 'a' and ' 6' res pectively. The a a and <7b are the RMS value for the timing 
residual R a (t) and Rb(t). One has (see lLee et al.l (|2008t ) for a similar calculation) 

m = — •£ (f° Km,. U ggL pA , (B2 , 

o-aCTfe x \io 647r/ g 7 + ck g • n a uj g + ck g ■ n b / 

where Tab is given by 

Vab = 1 - cos $ a - cos $ 6 + cos($ a - $ 6 ); 
$ a = — (Wg + cfc g • n a ); 

$ b = ^-(u g + ck g ■ h b ). (B3) 

Equation (|B2|) shows that the correlation function is composed of two major integrations, one over the GW frequency 
df g , and the one over solid angle dQ. For the massless graviton case, these two parts are not mixed, since fc g is linear 
in Lu g . But kg is non-linear in w g for the case of a massive graviton. The nonlinear dependence in Eq. ^ mixes the 
spatial integral with the frequency one. Thus the shape of C(9) depends on both the graviton mass and the frequency 
spectra of GWs, which also depends on the observation schedule. Thus, it is unlikely that one can integrate Eq. (IB2|) 
analytically. 



